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ABSTRACT 

We investigate the dynamics in the logarithmic galactic potential with an an- 
alytical approach. The phase-space structure of the real system is approximated 
with resonant detuned normal forms constructed with the method based on the 
Lie transform. Attention is focused on the properties of the axial periodic orbits 
and of low order 'boxlets' that play an important role in galactic models. Us- 
ing energy and ellipticity as parameters, we find analytical expressions of several 
useful indicators, such as stability-instability thresholds, bifurcations and phase- 
space fractions of some orbit families and compare them with numerical results 
available in the literature. 
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1. Introduction 

To determine salient features of the orbital structure of non-integrable potentials is 
a fundamental topic in galactic dynamics. Knowledge of existence and stability of periodic 
orbits of low commensurability is for example very important to clarify the issue of triaxiality 
in elliptical galaxies and a general understanding of the phase-space structure is very useful 
in applications of self-consistent models. The numerical approach is usually preferred due to 
the availability of rehable algorithms and powerful machines (Schwarzschild, 1979; Merritt & 
Valluri, 1996; Papaphilippou & Laskar, 1996 and 1998). However, in several circumstances, 
is useful to have some simple analytic clues concerning the relation between the form of the 
gravitational potential and the main families of orbits supported, as, for example, in the 
studies by Zhao et al. (1999). 

Although generically non-integrable, realistic galactic potentials show several properties 
of a regular behavior over wide energy ranges, with invariant surfaces ("tori") surrounding 
periodic orbits. Among these features, information about the stability properties of the 
main periodic orbits is of paramount importance, since the bulk of density distribution is 
shaped by the stars in regular phase-space domains around stable periodic orbits (Binney 
& Tremaine, 1987). In particular, for triaxial ellipsoids, periodic orbits along symmetry 
axes play a special role. An enormous effort has therefore been devoted to investigate 
families and bifurcations of periodic orbits, starting with the study of models based on 
perturbed oscillators (Contopoulos, 1970) and gradually exploring more realistic galactic 
potentials with numerical (Miralda-Escudc & Schwarzschild, 1989; Fridman & Merritt, 1997) 
and semi-analytical (de Zeeuw & Merritt, 1983; Scuflaire, 1995) approaches. A related 
important problem is that of torus construction in the regular domains (Binney & Spergel, 
1982; Gerhard & Saha, 1991; Kaasalainen & Binney, 1994a,b). 

Techniques based on the various versions of perturbation theory have been applied to 
several examples and with different degrees of approximation (for a review, see, e.g., Con- 
topoulos, 2002). One of the most powerful analytic tools is the normal form approximation of 
a non integrable system. Although the normal form approach is quite widespread in galactic 
dynamics, its use in studying stability of periodic orbits has not been as systematic as the 
theory could allow (Sanders & Verhulst, 1985). Aim of the present paper is to apply the Lie 
transform normalization method (Hori, 1966; Deprit, 1969; Dragt & Finn, 1976; Finn, 1984; 
Koseleff, 1994) to approximate the dynamics of the Binney logarithmic potential (Binney, 
1981). We compare the findings with that of Miralda-Escude & Schwarzschild (1989), who 
employ purely numerical techniques to implement the Floquet method and with that of Scu- 
flaire (1995), who studies the stability of axial orbits by solving the Hill-like perturbation 
equation with the Lindstedt-Poincare approach. Another example that is briefly mentioned 



-3- 



is provided by the galactic Schwarzschild (1979) potential with a comparison with the re- 
sults of de Zeeuw & Merritt (1983). These authors based their approach on the averaging 
procedure of normalization: it is therefore interesting a comparison with that method as 
well. 

To study the structure of phase space with a truncated normal form one may proceed 
in essentially two ways: the most general and exhaustive is that of determining the explicit 
form of periodic orbits and invariant tori that give the 'skeleton' of the phase-space of 
the system. A less general but easier approach in 2-degrees of freedom (DoF) systems, 
is that of determining fixed points and invariant curves on a surface of section. This is 
constructed with the aid of the approximate integral of motion provided by the normalization. 
Even in the case of linear orbital stability of the main periodic orbits, the first method is 
in general quite cumbersome and can be applied when the procedure of reduction to a 
single degree of freedom Hamiltonian system and the use of action-angle variables lead to 
a reasonably simple system of equations. This reduction is always possible with 2 DoF, 
both in the resonant and non-resonant cases (Gustavson, 1966), whereas it is in general no 
more possible in resonant 3-DoF systems. The second method is clearly less general but 
relies on straightforward geometric arguments related to the Hessian of a polynomial in its 
critical points and is, at least in principle, quite easy to implement. In this work we are 
going to apply both methods to perform the comparison mentioned in the paragraph above. 
There is another motivation to linger on both visions of the approximate dynamics: the 
normal form Hamiltonian is, in many respect, the "simplest" integrable relative of the non- 
integrable original one. However, this simplicity is attained at the price of a set of coordinate 
transformations that deform the phase-space variables; if one is interested in specific local 
properties or wants to make comparison with numerical simulations, the inversion to the 
original variables becomes mandatory. Therefore we want to deepen the relation between 
the 'clean' but deformed dynamics provided by the normalizing variables and the dynamics 
in the physical variables 'dirtied' by the inverse transformation. 

Even in realistic potentials, the periodic orbits along the axes of symmetry (axial orbits) 
are easily identified both as normal modes of the reduced system and as "central" fixed points 
on the surfaces of section. Therefore, we will limit the detailed evaluation of the stability 
characteristics in the parameter space to these axial orbits. However, both procedures we 
have followed are quite general and can be directly applied to all periodic orbits of sufficiently 
low commensurability. From the results obtained, we can state that the predictive power of 
the normal form ranges well outside the neighborhood in which the expansion of the original 
Hamiltonian is performed. It is rather related to the extent of the asymptotic convergence 
radius of the approximate integrals of motion, namely the extent at which the remainder 
of the series truncated at a given order is minimal (see, e.g., Efthymiopoulos, Giorgilli & 
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Contopoulos, 2004). The results suggest that there is still room for improvement up to a 
certain optimal order of truncation. Therefore, at least in the cases we have examined, the 
regular dynamics of a non-integrable system can be reproduced with very high accuracy. 
However, in concrete applications, the validity of the prediction has to be corroborated with 
an independent evaluation of the best suited resonant normal form for the problem at hand, 
always bearing in mind that, after the onset of chaos, these methods become ineffective. 

The plan of the paper is as follows: in Section 2 we recall the procedure of normalization 
as applied to reflection symmetry potentials and state a criterion for the choice of the best 
suited resonant normal form; in Sections 3 and 4 wc study the 1:1 and 1:2 resonances 
respectively; in Section 5 we compare our analytical results with those available in the 
literature; in Section 6 we sketch the conclusions. 

2. Analytical set-up 

The normal form approach to investigate Hamiltonian systems has a wide range of ap- 
plications in the astrophysical context. However, the method based on Lie transforms, which 
has several computational virtues with respect to the original methods (Gustavson, 1966; 
Contopoulos, 2002), has seldom been employed in galactic dynamics. Notable exceptions are 
the papers by Gerhard & Saha (1991) and Yanguas (2001) concerning perturbed isochrone 
models. 

In this section we recall the most relevant mathematical tools, with the aim of laying 
down the notations and the main formulas. For a detailed account of the method we refer 
to tetxbooks (Sanders & Verhulst, 1985; Meyer & Hall, 1992; Boccaletti & Pucacco, 1999). 

2.1. General 

Suppose the system under investigation is given by a Hamiltonian 

Hip,ci)^^ipl + pl) + Vix',y'), (1) 

with V a smooth potential with an absolute minimum in the origin and reflection symmetry 
with respect to both axes. We expand the potential in the form 

oo 

V{x,y;e) = Y,^'Vk{x,y) (2) 

k=0 
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and look for a new Hamiltonian given by 

oo 

K{P, Q;6) = J2 ^"^n(P, Q) = M;'H{p, q; e), (3) 



n=0 



where P, Q result from the canonical transformation 

(P,Q) = M,(p,q). (4) 

In these and subsequent formulas we adopt the convention of labeling the first term in the 
expansion with the index zero: in general, the 'zero order' terms are quadratic homogeneous 
polynomials and terms of order n are polynomials of degree n + 2. The linear differential 
operator Mg is defined by 

oo 
n=0 

The functions gn are the coefficients in the expansion of the generating function of the canon- 
ical transformation and the linear differential operator Lg is defined through the Poisson 
bracket 

L f ^ {g n + - - (6) 

^ ' dx dpx dy dpy dp^ dx dpy dy ' 

The exponentials in the definition of Mg are intended as the formal sum of a power series so 
that it gives rise to a near-identity coordinate transformation known as Lie series. Therefore, 
in practice, in the algorithm implemented to compute the Kn appearing in (I3l), what are 
really used are the differential operators 

T rn\ T m2 , , , T mn 

Mn= E (7) 

mi+2m2+...nmn=n 

By expanding also the right hand side of ([3]) in power series of e and equating the coefficients 
of the same order, one has the recursive set of linear partial differential equations 

Ko = Ho = 1{pI+pI) + Vo, 
K, = V, + M,Ho, 



Kn = Vn + MnHo + ^7=1 Mn-mVm , 
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It can be seen that each equation in the chain ([8]) depends only on quantities found in the 
previous hne. 'Solving' the equation at the n-th step consists of a twofold task: to find Kn 
and Qn- The unperturbed part of the Hamiltonian, Hq, determines the specific form of the 
transformation. In fact, the new Hamiltonian K is said to be in normal form if 

{Ho,K} = 0. (9) 

This condition is used at each step of the procedure to determine the functions gn in order 
to eliminate as much as possible terms in the new Hamiltonian. The only terms of which 
K is made of are those staying in the kernel of the operator Lh^ associated to Hq through 
definition ([6]). The procedure is stopped at some "optimal" order N and therefore in all 
ensuing discussion we refer to a "truncated" normal form. Hq must be considered a function 
of the new coordinates at each step in the process: according to ([9]), it is therefore an integral 
of the motion for the new Hamiltonian K. The function 

I = K -Ho (10) 

can be therefore used as a second integral of motion conveying approximate information on 
the dynamics of the original system. 

The normalizing transformation (jl]) leads to new coordinates which are continuous dif- 
ferentiable deformations of the original ones. For practical applications (for example to 
compare results with numerical computations) it is useful to express approximating func- 
tions in the original physical coordinates. Inverting the coordinate transformation, the new 
integral of motion can be expressed in terms of the original variables. Denoting it as the 
power series 

oo 

/ = 5^£"/., (11) 

n=0 

its terms can be recovered by means of 

n-l 

In = Vn-Kr,+ J2Mn-m[Vm-Im], n>l, (12) 

m=l 

that is obtained from (ITOl) by exploiting the nice properties of the Lie transform with respect 
to inversions (Boccaletti & Pucacco, 1999). 

It is important to remind in which respect the normal form K and the integral / provide 
approximations to the dynamics of the original system. The normal form is truncated at step 
A^: this means that in the new Hamiltonian a "remainder" of order 0(£:^'*'^) is neglected and 
the two functions K and I provide an (exactly) integrable system whose dynamics, in the 
new coordinates, is 0{e^) close to that of the original system. On the other hand, inverting 
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the canonical normalizing transformation, the power series f fTTj) is an approximate integral 
of motion of the original system (in the original coordinates) since 



That is, / fails to commute with the original Hamiltonian for terms of higher order. We 
see that in both cases the error made with the perturbative method is of the same amount 
since terms of the same order are neglected. However, it is not clear if, in the case of high 
orders and when a direct comparison is possible, both approaches lead exactly to the same 
predictions. Since in various applications one or the other approaches have pros and cons, 
one of the aim of the present work is just to examine and reconcile possible discrepancies. 

We remark that in all subsequent applications involving series expansions, the role of 
the perturbation parameter can also be played by the size of the neighbourhood of the origin 
where the Hamiltonian is considered. Therefore the powers of the parameter e are left in all 
expansion formulas just to indicate their order and are treated as unity in the computations. 



The model potential we will consider is the Binney logarithmic potential (Binney & 
Tremaine, 1987). The logarithmic potential has played a fundamental role in the description 
of galactic models displaying, with a very simple analytical expression, many realistic features 
of elliptical galaxies, in particular employing its singular scale-free form (Miralda-Escude & 
Schwarzschild, 1989) 



We cannot apply the standard normalization procedure as it is to a singular potential. In 
this way we are still not able to explore the phase-space structure of cuspy models that 
nowadays seem to be required by the observed properties of real galaxies (Merritt, 1999). 
However, selected classes of power law or other singular potentials can be put in a form 
suitable to the normalization algorithm by a proper coordinate transformation (Sridhar & 
Touma, 1997). Moreover, it is important to remark that the class of "weak" cusps (Dehnen, 
1993) displays many features of cored potentials as can be seen in the numerical exploration 
performed by Fridman & Merritt (1997). Therefore, the analysis presented here, even if 
restricted to potentials that can be Taylor expanded around a homogeneous density core, is 
a useful starting point. Enforcing the reflection symmetries of galactic potentials, we assume 
that each term in the expansion ([2]) can be written as a homogeneous polynomial of degree 
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2.2. Expanding "cored" galactic potentials 




(14) 
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2{k + l) of the form 



1 



V2k{x,y) 



2{k + l) 



j=0 



,2{k-j+l) 



A; = 0,1,2,... 



(15) 



A Taylor expansion has a convergence radius that can be finite and small if compared with 
the region of interesting dynamics. One could conservatively deduce that the normal form 
is a good appoximation of the dynamics only within the convergence radius of the potential. 
However, things are not so simple, because what really matters is the convergence of the 
normal form itself. In fact it may happen that a potential with infinite convergence radius, for 
example a perturbed oscillator without escape, has a normal form with a finite convergence 
radius: usually, this is related to the break up of regular dynamics and the transition to 
stochasticity. Clearly, in this case there is no hope to get reliable information on the dynamics 
beyond the stochasticity threshold from the normal form. On the opposite side, there is 
the possibility that the system, even if non integrable, displays regular dynamics almost 
everywhere and in this case the normal form can be used even outside the convergence 
radius of the original potential. The logarithmic potential offers the opportunity to explore 
all these issues (Kaasalainen & Binney, 1994b). 

The cored logarithmic potential can be written as 



The form written here is simplified by the choice of fixing the length scale (the "core radius" 
Rc) equal to one, but this is not a limitation due to the invariancc in both the length scale 
and the energy scale. With these units, the energy E may take any non negative value 



Lower values of q can in principle be considered but correspond to a non physical density 
distribution. Values greater than unity are included in the treatment by reversing the role of 
the coordinate axes, effectively extending the range until the upper limit 1/0.6 = 5/3. The 
series expansion of the logarithmic potential is 



y=^log(l + a;2 + ?/7g2). 



(16) 



< £; < oo. 



(17) 



The parameter giving the "ellipticity" of the figure ranges in the interval 



0.6 <q< 1. 



(18) 




(19) 
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so that the coefficients appearing in f[T5|l are 



'A; + 1\ (—1)''' 

a2U,k-j+i) = { J j q2ik-j+i) ■ (20) 



Those of lowest order are 



0(2,0) = = 1, «(0,2) = = V?^ (21) 

0(4,0) = -1, «(2,2) = -2/g^ 0(0,4) = -l/g^ (22) 
0(6,0) = 1, 0(4,2) = 0(2,4) = 3/5"^, 0(0,6) = (23) 

The natural setting in which one performs a low order normalization according to definition 
is therefore that of a perturbed quadratic Hamiltonian with a potential starting with the 
harmonic term and 'frequencies' given by (l2Ti) . For a generic value of q, the frequency ratio 

^ - . (24) 

is a real number: the frequencies are rationally independent and the terms in the normal 
form, that is those polynomials in the canonical variables commuting with 

^^0 = lipl +pI + + ivly'), (25) 

consist only of functions of the partial energies in the harmonic potential. It is customary 
to refer to the normal form constructed in this "Birkhoff" normal form (Birkhoff, 

1927). The presence of terms with small denominators in the expansion, forbids in general 
its convergence. It is therefore more effective to work since the start with a resonant normal 
form (Sanders & Verhulst, 1985), which is still non convergent, but has the advantage of 
avoiding the small divisors associated to a particular resonance. To catch the main features 
of the orbital structure, we therefore approximate the frequencies with a rational number 
plus a small "detuning" that we assume 0(5^) 

^ = g = ^ + e\ (26) 
UJ2 n 

We speak of a detuned {mm) resonance, with m+n the order oi the resonance. The algorithm 
to perform a resonant normalization generalizes that of the Birkhoff normalization in its 
ability to identify additional terms to be included in the normal form, since they cannot be 
eliminated with the sequence of canonical transformations. To this aim, the system must 
be in a form suitable to apply the resonant normalization procedure: we rescale variables 
according to 

X — > y/uhx, y — > ^/uT^V, Px — ^ -7^, Py — ^ -7^, (27) 
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in order to put the Hamiltonian in the form 

^ oo ^ fc+1 

n = -[{m + neH){pl + x') + n(pj + y')] + J] — ^ h,^,,,_,^,^x'^ y'^^-^^^^ (28) 
^ fc=0 ^^'^ + ^> i=o 

where we have used the same notation for the rescaled variables and 

na2(j-fc-j+i) 

Frequencies have been approximated as in fl26l) and the Hamiltonian is redefined according 
to the rescaling 

nH 

n:= = nqH. (30) 

As before, the new Hamiltonian K is said to be in normal form if 

{no,K} = (31) 

where, in agreement with 0301] . the unperturbed part of the Hamiltonian fl25p has been 
replaced by 

no = \[m{pl + x') + n{pl + y')] (32) 

and the procedure is now that of an ordinary resonant "Birkhoff-Gustavson" normalization 
(Gustavson, 1966; Moser, 1968) with two variants: the coordinate transformations are per- 
formed through the Lie series and the detuning quadratic term is treated as a term of higher 
order and put in the perturbation. 



2.3. Choice of the resonance 



The "skeleton" of the regular part of the phase-space of a non integrable system is 
framed by periodic orbits and invariant tori. In particular, knowing location and stability of 
periodic orbits allows one to gather a substantial piece of information concerning the whole 
dynamics. Moreover, determining the instability thresholds in terms of physical parameters 
(e.g., the energy), provides also clues on the extent at which stochasticity becomes important. 
The change in stability of a periodic orbit is connected to frequency ratios: commensurability 
of low orders between frequencies is the main trigger to stability-instability transition and 
the interaction of resonances provides chaos. In harmonic oscillators, frequency ratios are 
fixed: if their ratio is non-rational there is no resonance. For example, the quadratic part 
of the logarithmic potential is given by coefficients (|2T1) : with q in the range f|T8|) there is 
no reason to assume it to be a rational number. However, coupling between the degrees of 
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freedom due to the perturbation causes the frequency ratios to change. The system passes 
through resonances of order given by the integer ratios closest to the ratio of the unperturbed 
frequencies. Crossing the resonance, a stabihty change may occur and, in general, a new 
sequence of periodic orbits bifurcates (Merritt, 1999; Contopoulos, 2002). Given an arbitrary 
pair of unperturbed frequencies, it could seem better to approximate their ratio as close as 
possible with integers. However there is an argument on which a more effective choice can 
be based. A typical situation is that in which a family of periodic orbits becomes unstable 
when a low order resonance occurs between its fundamental frequency and that of a normal 
perturbation: the simplest case is given by an axial orbit that, depending on the specific form 
of the potential, can be unstable through bifurcation of loop orbits (1:1 resonance), "banana" 
orbits (1:2 resonance), "fish" orbits (2:3 resonance), etc. Therefore, a detuned low-order 
resonant normal form can be quite accurate in describing the corresponding bifurcations. 
The strategy can be either a systematic exploration of the hierarchy of resonances or a 
specific choice guided by some independent argument (energy range, ellipticity range, etc.). 
Clearly, since the frequency ratio (1261) is given by the parameter q that determines the shapes 
of the isopotentials, fixing a specific value of q in the range (fTSjl also indicates which resonance 
is expected to be closer. 

It must be emphasized that the structure of a resonant normal form is also affected by 
the symmetries of the original system. The normal form must preserve these symmetries and 
this in general also leads to a criterion for truncation. In the present instance of a double 
reflection symmetry, given a resonance ratio m/n, the normal form must contain at least 
terms of degree 2 x [m + n) (see, e.g. Tuwankotta & Verhulst, 2000). Therefore, the criterion 
we have adopted in this paper has been that of starting with the lowest order truncated 
normal form incorporating the symmetries of a typical galactic potential. A systematic 
investigation of the optimal order of truncation has recently been performed by Contopoulos 
et al. (2003) and Efthymiopoulos et al. (2004). Their results confirm the rapid decrease of the 
optimal order with the radius of the phase-space domain in which expansions are computed: 
we may conjecture that if we are interested in the global dynamics and accept a moderate 
level of accuracy, with this very conservative approach we can get reliable information up to 
the breakdown of the regular dynamics. Once particular orbits (or other interesting features) 
of the system are selected, more accurate details concerning the location, stability thresholds 
and so forth, can be derived with a higher order truncation. 
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3. 1:1 symmetric resonance and second order normalization 

We proceed by exploiting the normal form of the system with Hamiltonian fl28|l . po- 
tential f|T9l) and prescription fl29l) . starting with the simplest possibility: a Lie transform 
normalization truncated to the first non-null term in the normal form. Recalling the sym- 
metries of the potential represented by (I^Tj) . the first non trivial equation in the chain 
is 

K2 = V2 + M2H0, (33) 

so that we actually truncate at order 2. Applying standard methods (Meyer & Hall, 1992; 
Boccaletti & Pucacco, 1999) gives the following expression of the second-order normal form 
(Belmonte et al. 2006) 

= i5(Pl + X^) - I {qiPx' + Xr + \Py^ + 

--^ (Py\SPx^ + X^) + Y\Px^ + 3X2) ^ 4p^PyXY) , (34) 
lo 

Kt-^ = ^5(Pl + X') + h [q{Px' + Xr + l(Py2 + 

+k,{Px' + X'){Py' + Y'), (35) 

where the values of the coefficients in the expansion of the potential come from fl2Tl) and 
( I22I) and ki, k2 are rational numbers dependent on m and n. Eq.( l26|) has been used and the 
canonical variables P, Q are as in Eq.([3]). We see that this case behaves in the same way as 
in the first order averaging approach (de Zeeuw & Merritt, 1983; Verhulst, 1996): the 1:1 
resonance or all other resonances. This remark is another clue to the strategy delineated in 
subsection 2.3, that is to limit the global analysis to a normal form truncated to the first 
term incorporating a meaningful resonance. Therefore, for a preliminary investigation of 
periodic orbits with the lowest commensurability, that in the present instance is 

^ = 1 + 6^6, (36) 

U!2 

we can just exploit a 1 : 1 normal form truncated at 

In the following subsections we work out in detail this analysis to take it also as a record 
for the cases with higher order resonances. In fact, in these cases, the procedure goes along 
the same lines, but reporting computations in details is not so easy due to huge expressions. 
In Sect. 5.1 we will get more accurate numerical results going to order 6. A detailed study 
of the orbits in the 1:1 resonance case (including periodic and non periodic orbits) has been 
provided by Contopoulos (1965). 
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3.1. Orbit structure of the 1:1 resonance 

To investigate tlie phase-space structure of the system, it is convenient to write the 
normal form as 



K^^--^^ = J^ + J^ + £ 



- I (iJi + -J2 + %hJ2{2 + cos(20i - 262)) 
8 V q 3 



(37) 



where the action-angle variables are introduced according to 



X = y2Jicos^i, (38) 

= y2Jisin^i, (39) 

Y = y2Xcos^2, (40) 

Py = v/2^sin^2- (41) 

The structure of fl37|) displays the effect of the symmetries on the resonant part: angles appear 
only through the combination 26i — 202 and this shows why the symmetric 1:1 resonance 
can also be dubbed a "2:2" resonance. We now have an integrable system with dynamics 
generated by Hamiltonian (1371) with a second independent integral of motion 

no = S = Ji + J2. (42) 

This system is the simplest integrable approximation of the non-integrable dynamics in the 
logarithmic potential dominated by the lowest order resonance. In the normal form (l37|l . 
the perturbation parameter appears as a remind to denote terms of the same order and, 
according to the criterion stated above (see subsection 2.1), is put equal to unity in the 
computations. To get a quick overview of its structure we can use (1371) to identify the main 
periodic orbits. The procedure is the following (Sanders and Verhulst, 1985, sect. 7.4): we 
perform the canonical transformation to "adapted resonance coordinates" 

7/> = 2{0,-e2), (43) 

X = 2(6^ + 62), (44) 

Ji = i£ + n)/2, (45) 

J2 = {S-n)/2, (46) 

where £ is defined in f H2|) and 

n=Ji- J2. (47) 

Since x is cyclic and its conjugate action £ is the additional integral of motion, we introduce 
the "effective" Hamiltonian 

k = K^'-'\n,^-£,q), (48) 
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namely 

1 



K = 8 + ^{q - 1){S + TZ) + AiS"^ + n'^) + BSn + C{S^ - n^){2 + costfj), (49) 



with 



C = -i. (52) 

Considering the dynamics at a fixed value of S, we have that K defines a one-degree of 
freedom {ip,TZ) system with the following equations of motion 

tp = KTi = ^{q-l) + BS + 2{A-C{2 + cos^))n, (53) 
7^ = -K^ = c{S^- 7^2) sinV'. (54) 

Let us determine the fixed points of this system: these in turn give the periodic orbits of the 
original system. The right hand side of (15^ vanishes in one of the following three cases: 

I. For 7^ = ±£; 

II. For ^ = 0; 

III. For tfj = ±7r. 

In the first case, the right hand side of flS5]) vanishes when 

l-q-2[B±2{A-C{2 + cosip))]£ = (55) 
and the two periodic orbits 

n = £, Ja = 0, (Type la), (56) 

n = -£,Ji=0, (Type lb), (57) 

ensue. The orbit of Type la is the periodic orbit along the x-axis (long axial orbit) whereas 
the orbit of Type lb is the periodic orbit along the y-axis (short axial orbit). 

In the second case, the right hand of (153!) vanishes when 

q-l + 2B£ 8g-3(l + g)g 
^= A{3C-A) = 3{q-l) ' (^ = °) 
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where f l50] - [52l) have been used. This fixed point determines the "inchned" orbit 

J. = J2 = (Type II). (59) 

3(g- 1)' 3(g- 1) ' ^ ^ ^ ' 

Note that 

0<Ji,J2<^ (60) 
and this range determines the condition for existence of the orbit of Type II. 

In the third case, the right hand of ( !53l) vanishes when 



q-l + 2B£ 8g(g-l) + 3(l-g^)g 
^= A(C-A) = 3g^-2g + 3 ' = ^''^ 

The fixed point in fl6Tl) determines the elliptic ( "loop" ) orbit 

_ (3-^)£ + 4,fa^l) ^ ,(3,-l)g-4,fa^l) 

^ 3g2-2g + 3 ' ^ 3g2 - 2g + 3 ' ^ ' ^ ' 

The range (160!) still determines the condition for existence of the orbit of Type III. 



3.2. Stability analysis of the 1:1 resonance 

Let us now consider the question of the stability of the periodic orbits. As discussed 
in the previous section, the 1:1 resonant normal form essentially captures those features of 
the system characterized by the lowest order of commensurability between frequencies. In 
particular, it is able to describe the cases in which the nominal frequency characterizing a 
periodic orbit happens to become equal to that of a normal perturbation. The characteristic 
curve representing this equivalence in some suitable parameter space provides the stability- 
instability transition looked for. 

For orbits of Types II and III, an ordinary investigation of the equations of variations of 
the system ( I53f5^ is enough to perform the linear stability analysis (Contopoulos, 1978) in 
analogy with the Floquet method. However, in the case of axial orbits of Type I, action-angle 
variables have singularities on them and these affect also the adapted resonance coordinates. 
However, the remedy is quite straightforward: to use a mixed combination of action-angle 
variables on the orbit itself and Cartesian variables for the other degree of freedom. 

Let us start with the orbits of Type II and HI that, due to their non-singularity, are 
usually referred to as periodic orbits in general position. The system of differential equations 
for the perturbations {^ip, 671) is given by 
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The sign of the determinant of the Hessian matrix computed on the periodic orbit, PO, 

A = {k^^-KnTiK^^)\po (64) 

determines the fate of a perturbation: if A is negative it gives the frequency of bounded 
oscillating solutions 

Up = ^v^, (65) 

where the factor 1/g is due to the rescaling (l30l) . If A is positive it gives the characteristic 
exponent of the time evolution of a growing perturbation 

T = ±-Va. (66) 

q 

Therefore, with our 1:1 resonant normal form fH9|) . the condition for stability is 

A = AC^TZ^ sin^ ip - 2C{S^ + TZ^) cosip{A - C{2 + cosiIj))\po < 0. (67) 

Observing that each periodic orbit is identified by a fixed value of the pair {ip, TZ) and using 
the values of the constants in (ISUI) and (15^ . we see that the parameter space is spanned by the 
ellipticity q and the second integral of motion S: usually, characteristic curves in this space 
provide the instability threshold. We remark that the choice of the integral S as one of the 
coordinates in the parameter space is the simplest one and is usually adopted for a qualitative 
analysis (Sanders and Verhulst, 1985): however, for quantitative predictions and especially 
for comparisons with other approaches, a more natural choice would be the physical value of 
the energy E, of which S (but for a rescaling) is only a first order approximation. In general, 
using E in this framework is usually quite difficult for computational problems: this is one 
of the problems we mentioned above with exploiting the normal form in the normalizing 
coordinates. Using E with the original Hamiltonian and an approximate integral of motion 
is instead not only natural but also computationally easier. At the end of the subsection we 
will see how to overcome the problem in the particular case of orbits of Type I and how to 
express the characteristic curves in terms of E. 

An immediate application of condition fl^7|) concerns the inclined and loop orbits dis- 
cussed above. In the cases II and III, we get respectively 

A|^=o = 2{S^- 7^^) C{A - 3C) (Type II) (68) 

and 

A\^=^ = 2{S^ -n'^)C{C - A) (Type III). (69) 

From the obvious inequality £ > 71, A|^=o is always positive and A\^=t, is always negative: 
therefore the inclined (type II) orbit is unstable for every values of q and S, whereas the loop 
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(type III) orbit (when it exists, compare with conditions fl62|) and fl60|) ) is stable for every 
values of q and £. 

Our major concern is however to analyze the stability properties of Type I orbits, the 
axial orbits, also denoted as normal modes because each of them is identified by only one of 
the actions. Using action-angle variables on the normal mode and Cartesian variables on the 
normal bundle to it, the ensuing procedure is then first to determine the condition for the 
normal mode to be a critical curve of the Hamiltonian in these coordinates. Second, to assess 
its nature (Kummer, 1977; Contopoulos, 1978; Sanders and Verhulst, 1985, sect. 7. 4. 4): the 
condition is found by considering the function 

if ('^^ = K + /iHo, (70) 

where n has to be considered as a Lagrange multiplier to take into account that there is the 
constraint TCq = £. The Lagrange multiplier is found by imposing 

dif(^) = 0, (71) 

that is the total differential of (!70|) vanishes on the normal mode. Its nature is assessed by 
computing the matrix of second derivatives of if '^^^ : if the Hessian determinant of the second 
variation is positive definite the mode is elliptic stable; if it is negative definite the mode is 
hyperbohc unstable. In the case of the ?/-axis orbit of Eg. (1571) . good coordinates are X,Px 
and 

J = J2, (72) 
e = 02, (73) 

so that the periodic orbit is given by 

X = Px = 0, J = 8. (74) 
The terms in the normal form are then 

Ho = liX' + P'x) + J (75) 



(76) 



and 

I J [2{X^ + Pi) + (X2 - Pi) cos 2^ + 2XPx sin 2^] . 

It is straightforward to check that, in this case, Eq. flTTl) reduced to the periodic orbit defined 
by Eg. (1741) gives 

f. + l-- = 0, (77) 
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which allows us to find the required value of the Lagrange multiplier. With this result, the 
matrix of the second derivatives of K^^'^ on the periodic orbit is 



4(g- 1) -^[(2 + cos2^) - (3/g)] -Ssm29 

-£sm2e 4(g- 1) -^[(2-cos2^) - (3/g)] 



(7J 



The equation det[d^K^''\£)] = gives 

{q - l)(4g - 3S){Aq^ + 38- g(4 + ^)) = (79) 

and, according to the above recipe, for stability this polynomial in (g, £) must be positive. 
With q in the range 0181) as the independent parameter, the instability condition is 

'-^<S<^. (80) 

o — q 6 

It is interesting to remark that the instability threshold corresponding to the first inequality 
in (IHOl) coincides with the condition of existence of the loop orbit. In fact, from the obvious 
condition fl60|) and the values in (l62l) . we see that their validity is satisfied exactly by the first 
inequality in flHOl) : this is equivalent to state that at the point in the parameter space where 
the short axis orbits becomes unstable, the sequence of loop orbits bifurcates. The second 
inequality corresponds to a return to stability that actually disappears with the higher order 
treatment. 

Proceeding in the same way in the case of the x-axis orbit of Eq. (l56!) . analogous expres- 
sions can be obtained using non-singular coordinates J, 9, Y, Py and it is quite straightforward 
to check that, at this level of approximation, the Type la orbit is always stable. The pro- 
cedure above shows the amount of information that can be extracted from the 1:1 resonant 
normal form. In this respect, it is analogous to the first order averaging approach applied 
by de Zeeuw and Merritt (1983) in studying the Schwarzschild potential. 

Before going forward, we want to discuss how to give a more concrete meaning to the 
result enbodied by (IHOj) : this is possible if the physical energy E appears explicitly. According 
to the rescaling (!30l) . we assume that 

nE , , 

— (81) 

is the constant 'energy' value assumed by the truncated Hamiltonian K. In the present 
instance, n = 1 and 0^2 = 1/g, so that, on the ?/-axis orbit given by fl7^ . we have 

K = E- 1-^2 + _ ^ (82) 
8g 
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The dots are present to recall that a remainder has been neglected. The series fl82l) can be 
inverted to give 

g = q(E+^E' + ...) (83) 

and this can be used in the treatment of stability to replace S with E. With this substitution, 
the critical curve corresponding to the first inequality in (IHOi) is simply 

E{q)=2{l~q). (84) 

Expression (l82l) is also useful to find the oscillation frequency on the periodic orbit, 

which is needed if one is interested in computing either the frequency ratio 

r2 = ^, (86) 

where the perturbation frequency Up is defined in f l65l) . or the e-f aiding rate 

T2 = -, (87) 

K2 

where the time scale r is defined in fl66l). 



The analysis with a normal form truncated at higher order would provide more accurate 
predictions. However, for the time being, we want first to make the analogous analysis with 
the approximate integral, to better understand the relationship between the two approaches. 



3.3. Stability analysis with the approximate integral in the original variables 

Recalling the generic expression of the terms in the integral of motion (fT2l) . if we truncate 
to second order, we have 

I^'-'^ = Ho + e\V2-Ki''-'^). (88) 

This is the best approximate integral of motion of a symmetric perturbed 1:1 oscillator to 
order e"^ in the perturbation parameter: in fact, due to the symmetries of the problem, odd- 
degree terms are absent both in the normal form and in the approximate integral. A remark 
on how to interpret terms in (l88ll is necessary: even if not explicitly indicated, all of them 
must be intended as polynomials in the original variables {px,Py, x, y). In particular, K2^'^^ is 
the same as in fl34|) . where capital letter variables are simply replaced by the corresponding 
lower case letters. 
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To assess the stability of periodic orbits, we may proceed in the following way: starting 
for definiteness with the y-axis orbit, to use J^^-^) and the conserved energy E to construct 
an X —px Poincare section by means of the intersection of the function I^^''^\x, y,Px] E) with 
the y = hyperplane. The level curves of the function 

F = /(i^^)(a;,0,p,.;E) (89) 

allow us to determine the nature of critical points and invariant curves. Critical points {CP) 
on the surface correspond to periodic orbits in phase space (Contopoulos, 2002) and their 
nature gives the necessary information: they are either extrema (elliptic fixed point = stable 
periodic orbit) or saddles (hyperbolic fixed point = unstable periodic orbit) and this can be 
assessed by using the Hessian determinant 

{F,^v.Fxx-Fl^)\cp. (90) 

Clearly, for a periodic orbit, even the location of the critical points can be already quite 
difficult and this limits the generality of the approach. However, in the case of axial orbits, 
the approach is straightforward: the ?/-axis orbit, for example, coincides with the origin in 
this section. The second derivatives at the origin are 

4(g-l), (91) 
4(g-l) + 2E, (92) 
0. (93) 

The range of instability (namely, the range where the two non vanishing second derivatives 
have different sign) is 

E>2(l-g), (94) 

in agreement with (18^ . We remark that the results reported in Belmonte et al. (2006, 
cfr. Table 1), were obtained just by using both (ISO]) and (IMIl pointing out the discrepancy 
between the two predictions. An analogous procedure can be followed for the x-axis orbit by 
constructing a. y — py Poincare section and studying the level curves of the function 

F = I^'--'\0,y,Py-E) (95) 

obtained by means of the intersection of the function I^^'^\x, y,Py; E) with the x = hyper- 
plane. In agreement with the result of the previous subsection, the x-axis orbit is predicted 
to be stable using the present critical point analysis. 



FpxPx I CP 
Fxx\cp = 
FxpAcp = 
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4. 1:2 symmetric resonance and fourth order normalization 

The loop orbit bifurcates from the y-axis orbit at its instability threshold and is correctly 
described by the 1:1 (detuned) resonant normal form because its appearance coincides with 
the equality of the frequencies of the axial orbit and that of a normal perturbation. The 
'banana' orbits bifurcate from the x-axis when the frequency of the axial orbit falls to one 
half of that of a normal perturbation. This is the reason why the 1:1 resonant normal form 
does not detect any change of stability and we need the 1:2 (detuned) resonant normal form. 
A detailed study of the orbits in this resonance has been provided by Contopoulos (1963). 



4.1. Orbit structure of the 1:2 resonance 

In the case of the m = 1,?t, = 2 resonance in presence of reflection symmetries about 
both axes, we know that the normal form must be pushed at least to order 4, since it has 
to include terms of degree 2 x (m + n) = 6. Therefore, we have to perform a further step of 
normalization to include K4 in the normal form and the system to solve is now 

K2 = V2 + M^rCo, (96) 
K4 = Vi + M2V2 + M^Hq. (97) 

The expression of the normal form is quite involved (cfr. Belmonte et al. 2006), but we can 
exploit the change of variables to action-angle coordinates as deflned in fl38l) to see more 
easily its structure: 

= + 2J2 + e\25Ji - P2iJi, J2)) + e\P^{Ji, J2) + ^ cos(4^i - 2^2)), (98) 

8 

where the polynomials P2 and P3 are homogeneous of degree 2 and 3 respectively, 
P2 = I [qJ'i + ^J^^ + J1J2, 



and the frequency ratio (l26l) now is 

^ + sH. (99) 



0J2 2 

The adapted resonance coordinates are now deflned by 



i) = 461-262, 



(100) 
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X = 4^1 + 2^2, (101) 
Ji = £ + 2n, (102) 
J2 = 2£-n. (103) 

As before, knowledge of the orbit structure is obtained by investigating the fixed points of 
the one-dimensional dynamical system associated to the effective Hamiltonian 

K^'-^^ = K^'-^\n,ij-£,q), (104) 

where ^ 

7^ = -(2Jl- J2), (105) 

at fixed values of the integral of motion 

S = Uji + 2J2). (106) 



The study of existence and stability of periodic orbits proceeds in the same way as in 
the previous section. In the present case the normal modes (axial orbits) are given by 

n = 2£, J2 = 0, (Type la), (107) 
n = -£/2, Ji = 0, (Type lb). (108) 

The other periodic orbits are identified by the zeros of the system 

k!j^'-^^ = 0, (109) 
kI^'-^^ = 0. (110) 

In one case, (IllOp is satisfied hj ip = and the corresponding root of (1109^ 

7^ = nBi£, q), Ji = Ji(B), ■J2 = ■J2(B), (111) 

determines the "banana" orbit. In the other case, flllOp is satisfied hj ifj = n and the 
corresponding root of ( 1109^ 

n = TZa{£, q), Ji = Ji(A), J2 = J2(A), (112) 
determines the "antibanana" . From (11061) . their range is given by 

< JiiB,A) = ^T^bA^, q) + £<5£ (113) 

and 

< J2(B,A) = 2^ - ^b,a(^, q) < ^£. (114) 
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X 


= V2Jcos^ 


Px 


= \/2J sin 9, 


Y 


= Y, 


Py 


= V, 



4.2. Stability analysis with the normal form 

The procedure to determine the condition of stability of axial orbits leads to analyze 
critical curves of the modified function flTOj) where K is now given by fl98l) . From the analysis 
of the 1:1 resonance, we have seen that it cannot detect any instability threshold of the long 
axis periodic orbit and therefore we readily address this question in the framework of the 1:2 
resonance. Considering the x-axis (type la, Eq. fll07l) ) orbit, good coordinates are given by 

(115) 
(116) 
(117) 
(118) 

so that the periodic orbit is given by 

Y = V = 0, J = 58, (119) 

and 

= J + Y^ + V^. (120) 
Condition (I7T|) allows us to find the Lagrange multiplier 

= —(-32q + 2A£q - iOS^q + hlS'^q^). (121) 
16 

The equation det[i^'^^-*(£^)] = 0, obtained by computing the matrix of the second derivatives 
of K'^^'^ on the normal bundle to the periodic orbit flllOp . is 

(48-96g+24(3g-l)^+(26-153g+153g2)^2)(48_gg^+24(3g-l)£+(26-159g+153g2)£2) = 0. 

(122) 

whose zeros are 

^ ^ 3 - 9g ± v/3 V-23 + 187g - 432^^ + 306g3 

26 - 153g + 153g2 ' ^ > 



^A± = 4 — , . (124) 



3 - 9g ± V3 V-23 + 193g - 444g2 + 306g3 
26 - 159g + 153g2 

It happens that the conditions of existence in flll3|114l) are satisfied exactly when 

Sb+<S < Sb- (125) 

for the banana orbit and 

£a+<£ < £a- (126) 
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for the antibanana. We again see the relation between transition to instabihty of a normal 
mode and bifurcation of a resonant periodic orbit, since, outside the above ranges, the 
determinant A = det[i^^'^'*(£^)] is positive. As before, we want to make the analogous analysis 
with the approximate integral, to better grasp the relationship between the two approaches. 



4.3. Stability analysis with the approximate integral of motion 

We can develop the parallel approach of determining the nature of fixed points on the 
surface of section also for higher order approximate integrals. Using (fT2i) . if we truncate at 
order four, we have 

J(^^2) = Ho + e\V2 - K2) + e\V, - {g^, K^} - K,) (127) 

where 

g, = LplPyY-'^PlX--PxP^X-^PxX''-—P^Y-—PYX^Y--PxXY^-—PYY^ 
^ 24 ^ 16 ^ 3 ^ 16 32g ^ 24 6 32g 

(128) 

is the second order generating function determined in the first step of the normalization. This 
is the best approximate integral of motion of a symmetric perturbed 1:2 oscillator to order 
in the perturbation parameter. Let us again consider the x-axis orbit: the construction of 
the y — Py Poincare section and the study of the critical points of the function 

F = h{0,y,Py,E), (129) 

obtained by means of the intersection of the function I^{x,y,Py]E) with the x = hyper- 
plane, proceeds as above. Concerning the nature of the fixed point in the origin, we get 
quadratic inequalities in the energy analogous to those arising from (11221) : the discrepancy 
between the two methods is still of order 5^. The locus of points where the origin in the 
y — Py Poincare section changes from an extremum to a saddle are 

^ -8 + 10g2 ± y^228g - 887g2 + 830g3 + lOOg^ 
= ^ -76 + 165g^ ^ ^ 

^ -4 + 2g2 ± ^I80g - 659g2 + 614g3 + Aq^ 
= ^ 15(-4g + 7g2) " ^^^^^ 

These are clearly related to those found above but are not the same: to reconcile the two 
predictions, it is first necessary to give the relation between the new 'energy' £ and the 
physical energy E. But this is actually not enough: as seen in the 1:1 case, one must also 
perform a series expansion of the relations in the (g, E) plane and truncate it at the same 
order of the normal form. 
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According to the rescaling fl30l) . with n = 2 and uj2 = 1/q., on the x-axis orbit given by 
f fnOD . we have 

K = 2gJ - ^gj2 + - - 1)) + - = '^(lE- (132) 

The dots are present to recall that a remainder has been neglected. The series (11321) can be 
inverted to give 

j = ^£ = E+^-E'-^^(^^ + 17g(l - q)^ + ... (133) 

and this can be used in the treatment of stability to replace S with E. We limit ourselves 
to the lower bounds in the ranges in (11251) and (I126p . getting sufficient conditions for the 
bifurcation of, respectively, bananas and antibananas: 

i?B+=8(g-i)-y(g-i)^ (134) 
E^+ = 8(g-i) + y (g-i)^ (135) 

It can be checked by a long but easy computation that the above series coincide with those 
obtained by respectively expanding Eb+ in (I130p and Ea^ in (I13ip in q around q = 1/2, thus 
confirming the equivalence between the two methods. In the following section we present a 
more precise result obtained from a normal form truncated at higher order. 



5. Orbit structure of the logarithmic potential 

The theory discussed in the previous sections is applied to investigate the orbit structure 
of the logarithmic potential. In particular we are interested in simple recipes to predict the 
stability of the main periodic orbits. We can compare these analytical predictions with 
other results, mainly numerical, from the literature. In particular, we have chosen the work 
by Miralda-Escude & Schwarzschild (1989, MES in what follows) who made an accurate 
numerical exploration of the phase-space structure of the logarithmic potential: they give 
the existence parameter ranges of the main families of periodic orbits and the bifurcation 
ensuing from instability thresholds in two models (corresponding to ellipticity values q = 0.7 
and q = 0.9) determined by solving the perturbation equations with the Floquet method. 
The authors use the core radius Rc to parametrize the sequence of periodic orbits because 
they were interested in comparing the results with the case of the singular scale-free case 
Rc = 0: in fact they fix the energy E = in all their computations and vary Rc and q. 
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For our purposes, it is more natural to use the energy and the eUipticity as parameters. To 
compare our results with that in MES the conversion 

E^-logRc (136) 

must be used. We have numerically computed the instabihty thresholds of the axial orbits 
in the range 0.6 < g < 1, with results in agreement with MES where available. 

Another approach to the study of the stability of axial orbits in the logarithmic potential 
has been followed by Scuflaire (1995, SC in what follows) who solved the Hill-like equation 
of the normal perturbation to the periodic orbit by means of a Poincare-Lindstedt series 
expansion up to order 20. SC uses a, the amplitude of the axial orbit, as a parameter: the 
conversion to energy is given by 

E=^\og{l + a'), (137) 
E=Uog{l + (a/qy), (138) 

on the X-axis and y-axis orbits respectively. The approach in SC is analytical and perturbative 
but, being based on the theory of nonlinear differential equations with periodic coefficients 
(Jordan & Smith, 1999), is quite different from ours: it is therefore a very useful term of 
comparison because its results, even if restricted to axial orbits only, stem from a very high 
order perturbation analysis. 



5.1. Stability and bifurcation of closed orbits 

To get a higher precision than that obtained in the discussion of Sections 3 and 4, we have 
computed the 1:1 and 1:2 resonant normal form up to order 6 (degree 8 in the variables). A 
general analysis of these normal forms is quite cumbersome, but, concerning the axial orbits, 
they allows to get quite easily improved predictions of the stability thresholds. 

In the case of the y-a.xis orbit, using the same coordinate as in Subsection 3.2, the 
normal form on the periodic orbit ( Ji — 0, J2 = E) is 

%q 192g2 I024g3 ^ ^ ' 

As a consequence, the upgraded evaluation of the frequency on the ?/-axis orbit is 

IdK If 3 ^ 29 . 55 , , , 
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These series generalize f l82|) and (l85l) respectively. As before, we actually prefer to have 
expressions in terms of the 'real' energy E, so we invert fll39p to get 

S = qE(l + -E + —E' + -^E^ + ..) (141) 
^ V 8 192 1024 J ^ ' 

and substitute in (11401) . obtaining 

Computing the determinant of the matrix of the second derivatives of the modified function 
flTOj) . det[(fK^^'\£)] and putting it equal to zero, gives 

(q - l)(64g2 - A8q£ + (29 - 17q)£^){6Aq'^{l + q) + 16g(3 - q)£ - (g^ -2q + 29)£^) = (143) 

generalizing (17^ . Substituting (I14ip and solving for E we get 

^;(g)=2(l-g) + (l-g)2-^(l-g)^ (144) 

obtaining a prediction of the stability-instability transition of the short axis orbit (with 
subsequent bifurcation of the loop orbit) up to third order in the (g, E) plane. The values 
obtained from this relation are shown in Fig.([T]), continuous line, and are compared with the 
numerical results with the Floquet theory (dotted line) showing a very good agreement: the 
values obtained by MES (some of the dots in the figure) are 0.21 in the case q = 0.9 and 
0.72 in the case g = 0.7. Moreover, the expansions (11401) and (11441) coincide up to the same 
order with those reported in SC when the conversion from amplitude to energy is performed 
using ffM . 

We may also compute the ratio fl86|) between the frequency of the perturbation and that 
of the periodic orbit, where the perturbation frequency is associated to the determinant in 
the left side of (I143p like in fl^ and the e- folding rate as defined in fl57|) : as expected, at low 
energy, the ratio has limit q, since cjp ^ 1, ^2 — > 1/g, whereas, at the transition, the ratio 
coincides with the resonance value 1. In Fig. ([2]) the two quantities are plotted for q = 0.9, in 
Fig. ([3]) they are plotted for q = 0.7. They can be compared with the corresponding figures 
in MES (bottom panels in their figg.7 and 8) where they have been obtained numerically. 
The return to stability at higher energies mentioned in Sect. 3 (cfr. the second inequality 
in (I80p ) does not occur in the logarithmic potential. This incompleteness appears also in 
the treatment made by de Zeeuw & Merritt (1983) of the Schwarzschild potential using first 
order averaging. This problem is solved by our higher order approach, because the other 
roots of fll43p are either negative or complex. 
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In the case of the x-axis orbit, we hmit ourselves to provide results without derivation. 
The 1:2 normal form on the periodic orbit (Eq. fl98l) with J2 = 0) is 

As a consequence, the upgraded evaluation of the expansion of the action in terms of the 
true energy is 

The frequency on the x-axis orbit computed through 

..(JO-^f, (147) 
and expressed in term of the upgraded expansion (11461) , is 

K, = l--E + —E^ + —E\ (148) 
4 64 256 ^ ' 

Computing the determinant of the matrix of the second derivatives of the modified function 
(ITOjl on the periodic orbit, A = det[(PK^^^^] and putting it equal to zero, gives, in this case, 
an equation of sixth degree: the two roots in addition to those in (11301) and (11311) . in analogy 
to the roots Eb- and Ea-, does not put further constraints on the conditions sufficient for 
instability of the axial orbit. Therefore, we use the upgraded version of the relevant roots 
and express them as series in powers of g — 5, to get 

Eb^ = 8 (g -I) - f {q -\f + ^ (g -\f , (149) 
Ea+ = 8 (g -1) + I (g -\f + ^ (g -if , (150) 

obtaining a prediction of the stability-instability transition of the long axis orbit (with 
subsequent bifurcation of the banana orbit) up to third order in the (g, E) plane. The values 
obtained from these relations are shown in Fig.(jl]): the continuous and dotted lines are 
respectively the relations in (11491) and (11501) . and are compared with the numerical results 
with the Floquet theory, namely bifurcation of the banana (dash-dotted line) and bifurcation 
of antibanana (dashed line). The values obtained by MES are respectively E = 1.52 and 
E = 4.29 in the case g = 0.7, whereas, in the case g = 0.9, only the first transition at 
E = 3.62 is available. Moreover, the expansions (I148H1501) coincide up to the same order 
with those reported in SC when the conversion from amplitude to energy is performed using 
(II37D- 
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We may also compute the ratio 

ri = ^, (151) 

between the frequency of the perturbation and that of the periodic orbit, where the pertur- 
bation frequency is associated to the determinant in the left side of (11221) like in (!65|) and 
the e- folding rate as defined in (!87|) . but with ki in place of K2- in this case, at low energy, 
the ratio has limit 1/q, since tUp — > 1/g, ki — > 1, whereas, at the transition, the ratio has the 
resonance value 2. In Fig.([S]) the two quantities are plotted for q = 0.9, in Fig.([S]) they are 
plotted for q = 0.7. They can be compared with the corresponding figures in MES (upper 
panels in their figg.7 and 8) where they have been obtained numerically. 

A demanding test of an analytical approach like that based on normal forms is that 
of investigating higher-order boxlets: actually the procedure is analogous to that above, 
but the number of terms in the normal form can be very high. We have limited ourselves to 
investigating the bifurcation of fish orbits. To this aim we have to compute the 2:3 symmetric 
normal form which must include terms of degree at least 2 x (2 + 3) = 10. Once obtained 
we can express it in adapted resonance coordinates 

^ = 6^1 - 4^2, (152) 

X = 6^1+4^2, (153) 

Ji = 2£ + 3n, (154) 

J2 = 3£~2n. (155) 

and, as before, knowledge of the orbit structure is obtained by investigating the fixed points 
of the one- dimensional dynamical system associated to the effective Hamiltonian 

K^'-^^ = K^^--'\n,t,S,q). (156) 

The periodic orbits are identified by the zeros of the system 

k!^'-^^ = 0, (157) 
Kf-^^ = 0. (158) 

In particular, fll58p is satisfied by = and the corresponding root of (11571) 

n = nF{£,q) (159) 

determines the "fish" orbit. Using this relation and fll54p in the existence condition < 
JiF < 2£^, gives an equation whose roots are related to the bifurcation of the fish orbit from 
the normal mode. Using as usual the true energy in place of the fictitious £, we get the 
fourth degree equation 



Y,c^{q)E' = 0, (160) 



4=0 
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with coefficients 

Co = -6 + 9q, 

Cl = O — , 



C2 = 
C3 = 
C4 = 



4 

381 54639(7 677727g2 14337g3 12393?'* 



500 4000 16000 256 512 ' 

637 3186g 705807g2 131787^3 663399-* 
250 ~ 125 8000 1024 ^ 1024 ' 

18689 4279197g 28267497?^ 15582699q3 1062755397^4 1370678193g5 6134535^6 5688387g^ 



12800 256000 7168000 102400 3276800 4096000 32768 131072 

Equation (11601) has a complex pair and a real pair of solutions: the smaller positive root can 
be used for our purpose as a prediction of the bifurcation energy. For q = 0.7 and q = 0.9, 
we get Ep = 0.2 and Ep = 1.7 respectively. The bifurcation energies numerically found 
in MES are Ep = 0.21 and Ep = 2.28 respectively: clearly, the disagreement between the 
theoretical and numerical results rapidly grows with detuning. 



5.2. Phase space fraction occupied by boxlets 

An analytical estimate of the fraction of phase space occupied by each orbit family can 
be made with the tools constructed so far. The method is based on the computation of the 
area on the surface of section associated to the given orbit family. In this respect the estimate 
is affected by an error due to neglecting a factor depending on the period of individual orbits: 
this factor in general doesn't allow us to assume simple proportionality between volumes in 
phase space and areas on the surface of section (Binney, Gerhard & Hut, 1985). However, 
in our case, the periods of the families at hand belong to a range small enough to produce 
a negligible effect. To illustrate the method in its simplest form, we use the lowest order 
theory applied to loop orbits in Sect. 3 and to bananas in Sect. 4. 

For the loop orbits we have that fl62l) give the values of the actions in terms of S and q. 
As usual, to compare with quantities of the 'real' system, we need to convert to the physical 
energy and this can be done by means of the K on the loop 

3q^ — 2q + 3 

Inverting and expanding, this gives the integral S = S{E) in terms of the true energy. An 
important point to remark on, is that the J's appearing in the general expression of the 
normal form, e.g. K^^'^^ of fl37|) . are not true actions along a given family. However, the J's 
in (1621) are true actions: in particular, at the bifurcation of the loop from the axial orbit, 
they are the actions of the orbits asymptotic {homoclinic) to the y-axis orbit, the "ffist" loop 
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and the "last" box respectively. Therefore, Ji(£^)|7// is the area inside the homoclinic orbit 
on the X — surface of section and J2(-E)|7// is the area inside the homoclinic orbit on the 
y — Py surface of section. The ratios 

r _ Ji{.E)\iu , _ J2{.E)\iii , . 

jLoop c/T-iN ' JBox c f 1 \^^'^) 

can therefore be used as an approximation of the fraction of phase space occupied by the loops 
and the boxes respectively. Substituting into (162!) and dividing by S{E), their expressions 
are 

_ 2(-3 + 3g - + 5q^) + E{9 - 9q + iW - 3g') , . 

■^"-""^ (3-2g + 3g2)(-2(l -g)2 + E(3-2g + 3g2)) ^ ' 

and 

_ g(10 - lOg + - Qq^ + E(-3 + llg - + V)) 
}box- ^3_2g + 3g2)(_2(l -g)2 + E(3-2g + 3g2)) " 

In Figg. (!7f8l) the two quantities are plotted for q = 0.7 and q = 0.9 respectively and a 
comparison can be made with the numerical estimate of fioop in MES (cfr. their Fig. 10) 
computed by applying an approximate version of the recipe by Binney et al. (1985). We 
remark that with their approach, MES are not able to compute the fraction pertaining to 
boxes, since their method does not use families associated to normal modes. 

Concerning banana orbits, we may follow an analogous line of reasoning, obtaining 

J Ban — ^^j^^ , UD£)j 

where Ji(b), introduced in (lllip . can be interpreted as the area inside the orbit homoclinic to 
the unstable fixed point on the x—px surface of section. In Fig. ([7]) this quantity is plotted for 
q = 0.7 whereas, for q = 0.9, its application is meaningless because it concerns exceedingly 
high energies. Comparing with the numerical estimate of fsan in MES (cfr. the left panel of 
their Fig. 10), we have to take into account that our analytical estimates originate from two 
independent systems, the two normal forms based on the 1:1 and 1:2 resonances. Therefore, 
we cannot expect that the occupancies of each family sum up to give the whole phase space. 
A suitable normalization would be necessary to adjust the proportions and could account 
for the numerical result. However, since the normalization process contains a certain degree 
of arbitrariness, we let the quantity as they are and use them only as relative measures. We 
come again to this point in the next subsection. 
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5.3. Surface of section for singular logarithmic potential 

It is tempting to try to extract information concerning tlie scale-free singular limit of tlie 
logarithmic potential from our analytical setting based on series expansions. Formally, this 
operation should be hindered by the lack of a series representation of the singular logarithmic 
potential. However, we may nonetheless 'force' our approximate integrals of motion to play 
their role in the singular limit too. If we try our chance by constructing e.g. the y — Py 
surface of section with the same procedure as in Subsect. 4.3, using now 

lipl +pI) + 1 ^og{x' + y'/q') = (166) 

to eliminate p^, we get, quite surprisingly, acceptable results. In Fig. ([9]) we see the section 
obtained by using the approximate integral /(^-^^ of flHHl) for q = 0.7: it gives the family of 
loops around the stable periodic orbit at y ~ 0.56 and can be compared with Fig.l in MES. 
In Fig. flTUl) we see the section obtained by using the approximate integral /(^-^^ of fll27p . 
again for q = 0.7: it gives the family around the stable banana at ?/ ~ 0.16 and boxes around 
it. A comparison with the same Fig.l in MES shows that also the banana is located quite 
well on the section. However, from these considerations emerge the main limitation of the 
approach based on resonant normal forms: each resonance captures only one feature of the 
system under study. For example, comparing our Figg. ipiTOl) with a numerically computed 
section like that in Fig.l in MES, we see that, in the 'true' system, the island around the 
banana 'eats' a substantial part of the island around the closed loop. A reliable description of 
this phenomenon is not possible in the framework of a single-resonance theory. This remark 
is evidently related to the interpretation of the previous results concerning the computation 
of the fractions in phace space. 



6. Conclusions 

We have applied the method of resonant detuned normal forms to investigate relevant 
features of a non-integrable potential of interest in galactic dynamics. The analytical set up 
consists in constructing a new integrable system (the 'normal form') by means of a sequence 
of canonical transformations with the method of the Lie transform. The algorithm can be 
put in an efficient and powerful form (Boccaletti & Pucacco, 1999; Giorgilli, 2002) which can 
handle quite high-order expansions. We have shown how to exploit resonant normal forms to 
extract information on several aspects of the dynamics of the original system. In particular, 
using energy and ellipticity as parameters, we have computed the instability thresholds of 
axial orbits, bifurcation values of low-order boxlets and phase-space fractions pertaining to 
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the families around them. We have also shown how to infer something about the singular 
limit of the potential. 

As in any analytical approach, this method has the virtue of embodying in (more or less) 
compact formulas simple rules to compute specific properties, giving a general overview of the 
behavior of the system. In the case in which a non-integrable system has a regular behavior 
in a large portion of its phase space, a very conservative strategy like the one adopted in this 
work provides sufficient qualitative and quantitative agreement with other more accurate but 
less general approaches. In our view, the most relevant limitation of this approach, common 
to all perturbation methods, comes from the intrinsic structure of the single-resonance normal 
form. The usual feeling about the problems posed by non-integrable dynamics is in general 
grounded on trying to cope with the interaction of (several) resonances. Each normal form is 
instead able to correctly describe only one resonance at the time. However, we remark that 
the regular dynamics of a non-integrable system can be imagined as a superposition of very 
weakly interacting resonances. If we are not interested in the thin stochastic layers in the 
regular regime, each portion of phase space associated with a given resonance has a fairly 
good alias in the corresponding normal form. An important subject of investigation would 
therefore be that of including weak interactions in a sort of higher order perturbation theory. 
For the time being, there are two natural lines of developement of this work: 1. to extend 
the analysis to cuspy potentials and/or central 'black holes'; 2. to apply this normalization 
algorithm to three degrees of freedom systems. 
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Fig. 1. — Stability thresholds of the short y-axis orbit: analytical (Eq. 11441 continuous line); 
numerical solution with the Floquet method (dotted line). 
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Fig. 2. — Frequency ratio and e-folding rate for the short y-axis orbit at ^ = 0.9. 



-38- 



1 



0.8 - 



o . 6 - 



''-10.4- 



0.2 - 




0.25 0.5 0.75 1 1.25 1.5 



Energy 

Fig. 3. — Frequency ratio and e-folding rate for the short y-axis orbit at g = 0.7. 
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Fig. 4. — Stability thresholds of the long x-axis orbit: analytical (banana bifurcation, Eq. 
11491 continuous line and antibanana bifurcation, Eq. 11501 dotted line) and numerical (dash- 
dotted and dashed line respectively). 
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Fig. 5. — Frequency ratio and e-folding rate for the long x-axis orbit at q — 0.9. 
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Fig. 6. — Frequency ratio and e-folding rate for the long x-axis orbit at ^ = 0.7. 
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Fig. 7. — Fraction of phase space occupied by low order boxlets a.t q = 0.7: loops (continuous 
line); boxes (long dashed line); bananas (short dashed line). 
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Fig. 8. — Fraction of phase space occupied by low order boxlets at g = 0.9: loops (continuous 
line); boxes (dashed line). 
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Fig. 9. — Zero-energy y — Py surface of section of the singular logarithmic potential with 
q — 0.7, computed with the 1:1 resonant normal form. 
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Fig. 10. — Zero-energy y — Py surface of section of the singular logarithmic potential with 
q — 0.7, computed with the 1:2 resonant normal form. 



